A high resolution foreground cleaned CMB map from WMAP 
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We perform an independent foreground analysis of the WMAP maps to produce a cleaned CMB 
map (available online) useful for cross-correlation with, e.g., galaxy and X-ray maps. We use a 
variant of the Tegmark & Efstathiou (1996) technique that assumes that the CMB has a blackbody 
spectrum, but is otherwise completely blind, making no assumptions about the CMB power spec- 
trum, the foregrounds, WMAP detector noise or external templates. Compared with the foreground- 
cleaned internal linear combination map produced by the WMAP team, our map has the advantage 
of containing less non-CMB power (from foregrounds and detector noise) outside the Galactic plane. 
The difference is most important on the the angular scale of the first acoustic peak and below, since 
our cleaned map is at the highest (12.6') rather than lowest (49.2') WMAP resolution. We also 
produce a Wiener filtered CMB map, representing our best guess as to what the CMB sky actually 
looks like, as well as CMB-free maps at the five WMAP frequencies useful for foreground studies. 

We argue that our CMB map is clean enough that the lowest multipoles can be measured without 
any galaxy cut, and obtain a quadrupole value that is slightly less low than that from the cut-sky 
WMAP team analysis. This can be understood from a map of the CMB quadrupole, which shows 
much of its power falling within the Galaxy cut region, seemingly coincidentally. Intriguingly, both 
the quadrupole and the octopole are seen to have power suppressed along a particular spatial axis, 
which lines up between the two, roughly towards (l,b) ~ ( — 110°, 60°) in Virgo. 



I. INTRODUCTION 

The release of the first results [1-18] from the Wilkin- 
son Microwave Anisotropy Probe (WMAP) constituted 
a major milestone in cosmology, laying a solid founda- 
tion upon which to found the cosmological quest in com- 
ing years. Although much of the attention in the wake 
of the WMAP release has focused on "sexy" issues like 
the power spectrum and its cosmological implications, 
the primary stated science goal of WMAP is to produce 
maps. Indeed, one of the qualitatively new types of re- 
search made possible by WMAP involves taking advan- 
tage of this spatial information by cross-correlating the 
maps with other cosmological templates such as galaxy 
[19], x-ray [20, 21], infrared and lensing maps [22], which 
can reveal interesting signals ranging from the Late In- 
tegrated Sachs- Wolfe effect to the SZ effect and lensing 
[23, 24]. 

Many such future studies will be looking for signals of 
modest statistical significance, so it is important to quan- 
tify and minimize unwanted signals in the map due to 
foreground contamination and detector noise. Accurately 
understanding the foreground signal is also important for 
the interpretation of the WMAP early reionization detec- 
tion [1, 15, 17, 25-27] and for the interpretation of the low 
WMAP quadrupole [1, 6, 17], since Galactic foregrounds 
are most important on large angular scales [28, 30]. The 
WMAP team has already performed a careful foreground 
analysis [3] combining the five frequency bands into a 
single foreground-cleaned map, shown in Figure 1 (top). 
Given the huge effort that has gone into creating the 



spectacular multifrequency maps, it is clearly worthwhile 
to subject them to an independent foreground analysis. 
This is the purpose of the present paper, which we will 
argue not only corroborates the findings of the WMAP 
team, but also makes some further improvements that we 
believe are useful. 

The main goal of this paper is to remove foregrounds, 
not to understand or model them. For reviews of fore- 
ground modeling issues, see, e.g., [3, 28-41] and refer- 
ences therein. There is a rich literature of techniques 
for foreground removal. The work most closely related 
to the present paper is that done in preparation for 
the Planck mission [33-35], developing multipole-based 
cleaning techniques and testing them on simulations. 



II. ANALYSIS AND RESULTS 

A. The problem 

A key goal of the CMB community is to measure the 
function x(r), the true CMB sky temperature in the sky 
direction given by the unit vector r. The WMAP team 
have observed the sky in five frequency bands denoted 
K, Ka, Q, V and W, centered on the frequencies of 22.8, 
33.0, 40.7, 60.8 and 93.5 GHz, respectively, producing five 
corresponding maps (Figure 2, left) that we will refer to 
as yi, i — 1,...,5. In practice, each of these maps are 
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FIG. 1. The linearly cleaned WMAP team map (top), our Wiener filtered map (middle) and our raw map (bottom). Here and 
throughout, all maps are shown in Mollweide projection in Galactic coordinates with the Galactic center (I, b) = (0, 0) in the middle and 
Galactic longitude I increasing to the left. These low-resolution images do not do justice to the WMAP data or our method; full res 
versions are available at http : / /www. hep.upenn.edu/ max/ wmap.html. 
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discretized into n = 12 x 512 2 = 3,145,728 HEALPix 1 
pixels [42, 43], so y, is an n-dimcnsional vector. However, 
since the maps are more than adequately oversamplcd 
relative to the beam resolution, it is equivalent and often 
simpler to think of them simply as five smooth functions 
yi (?) . Conversely, we will occasionally write the true sky 
as a pixelized vector x. 

These maps are related to the true sky x by the affinc 
relation 

y^AjX+rii, (1) 

where the matrix A, encodes the effect of beam smooth- 
ing and the additive term rij is the contribution from de- 
tector noise and foreground contamination, collectively 
referred to as "junk" below since it complicates the re- 
covery of x. Defining the larger matrices and vectors 

/Ai\ /yi\ /m\ 

A = : , y= : , n= : , (2) 
VA 5 ; \y 5 J \n 5 J 

we can rewrite the system of equations given by (1) as 

y = Ax + n, (3) 

a set of linear equations that would be highly over- 
determined by a factor five if it were not for the presence 
of unknown junk n. 

Foreground removal involves inverting this overdeter- 
mined system of noisy linear equations. The most general 
linear 2 estimate of x of the true sky x can be written 

x = Wy (4) 

for some n x (5n) matrix W. We will require that the 
inversion leaves the true sky unaffected, i.e., that the 
expected measurement error (x) — x is independent of 
x. Bennett et al. [3] refer to this property as the in- 
version having unit response to the CMB. Methods in- 
volving maximum-entropy reconstruction or some form 
of smoothing typically lack this property. Substituting 
equation (3) into equation (4) shows that this require- 
ment corresponds to the constraint 

WA = I. (5) 

B. The mathematically optimal solution 

Which choice of W gives the smallest rms errors from 
foregrounds and detector noise combined? Physically dif- 
ferent but mathematically identical problems were solved 



1 The HEALPix package is available from http://www.eso.org/ 
science/healpix/. 

2 In addition to simplicity and transparency, linear methods have 
the advantage that the noise properties of x can be readily cal- 
culated from those of the input maps. 



in a CMB context by [44, 45] , showing that the minimum- 
variance choice is 

x = [A*N _1 A] _1 A*N _1 y! (6) 

where N = (nn*). For an extensive discussion of dif- 
ferent methods proposed in the literature and the rela- 
tions between them, see [29]. Although optimal, this 
method is unfortunately unfeasible for the WMAP case, 
for two reasons. First, it requires the inversion of the 
(5n) x (5n) matrix N. Although the detector noise con- 
tribution to this matrix is close to diagonal for WMAP, 
the foreground contribution is certainly not. Second, it 
requires knowing the matrix N, which has many more 
components than there are pixels in the map. 



C. The WMAP team solution 

In producing their internal linear combination (LLC) 
map, the WMAP team adopt a simpler approach [3], first 
smoothing all five maps to a common angular resolution 
of 1° and then performing the cleaning separately for each 
pixel. The smoothing implies that Aj = I (if we redefine 
x to be the true sky smoothed to 1°), and equation (4) 
reduces to the simple form 

x(f) = ^ Wj y 2 (r) (7) 

for five weights Wi that according to equation (5) must 
sum to unity. The WMAP team chose the weights that 
minimize the rms fluctuations in the cleaned map x, using 
12 separate weight vectors for 12 disjoint sky regions. 

Although this method works well, it can be improved 
by allowing the weights to depend on angular scale (i.e. 
on harmonic number t) as well as on Galactic latitude. 
This has two advantages. First, the angular resolution 
is limited not by that of the lowest resolution channel 
(the K-band FWHM is 49.2'), but rather by that of the 
highest resolution channel (the W-band FWHM is 12.6'). 
Second, as shown by [28], letting the weights depend on 
angular scale produces a cleaner map by taking into ac- 
count that the frequency dependence of the unwanted 
signals varies with scale: at large scales, galactic fore- 
grounds dominate, whereas at small scales, detector noise 
dominates. 



D. Our solution 

In this paper, we will take an approach intermediate 
between the two described above, aiming to approxi- 
mate the optimal method while staying within the realm 
of numerical feasibility Our method is essentially that 
of Tegmark & Efstathiou [28], implemented to make 
it blind and free of assumptions both about the CMB 
power spectrum and about foreground and noise prop- 
erties. The only assumption, which is crucial, is that 



■5 




FIG. 3: The weights W( used to create the internal lin- 
ear combination map of the WMAP team are independent 
of angular scale. The figure shows the weights w e = 
(0.109, -0.684, -0.096, 1.921, -0.250) used outside of the galactic 
plane. 




FIG. 4: If there were no foregrounds and equal noise in the five 
input maps, then equal weighting at low £ would give way to fa- 
voring the highest resolution bands at high £. This example uses 
the forecast WMAP beam and noise specifications from [29] rather 
than the actual ones. 



the CMB has the Blackbody spectrum determined by 
the COBE/FIRAS experiment [60]. The method strictly 
respects the requirement of equation (5), which is most 
easily seen by verifying that each of the steps described 
below do so individually. The gist of our method is to 
combine the five WMAP bands with weights that depend 
both on angular scale and on distance to the Galactic 
plane (we subdivide the sky into 9 regions of decreasing 
overall cleanliness). We begin by describing the angular 
scale separation, then turn to the spatial subdivision in 
Section II F. 

We perform our cleaning in multipole space as in [28], 
and go back and forth between spherical harmonic coef- 
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FIG. 5: The optimal WMAP weights forecast by [29] for the 
middle-of-the-road foreground model from [29]. 




FIG. 6: The actual weights we use for the 3rd cleanest of the 9 
sky regions shown in Figure 8. The dirtier the sky region, the more 
aggressive the weighting becomes, using large negative and positive 
values to subtract foregrounds. 

ficients 

4m = I r, m (r)*x 4 (r)dO (8) 
and real-space maps 

x i^) = J2 a im Y em{r) (9) 

using the HEALPix package [42, 43] with £ max = 1024. 
Since this employs a spherical harmonic transform al- 
gorithm using fast Fourier transforms in the azimuthal 
direction [46], each transformation takes only about a 
minute on a Linux workstation. Each cleaned map is 
defined by 

<^ = X>1% (10) 
«=i a t 
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for some set of five-dimensional weight vectors Wf, where 
B\ is the beam function for the i th channel from [14]. 
(There are 4 W-band maps, 2 V-band maps and 2 Q- 
band maps; we combine these into single maps at each 
frequency by straight averaging and therefore average the 
corresponding beam functions as well.) When computing 
our final cleaned map in real space, we multiply a\ m by 
.Bf in equation (9) so that it has the beam corresponding 
to the highest-resolution WMAP band. 

To gain intuition for the weight vectors wi that specify 
a cleaned map, we have plotted them for four interesting 
cases in figures 3, 4, 5 and 6. Figure 3 corresponds to the 
weighting used by the WMAP team for the region away 
from the galactic plane, and is simply independent of 
I. To recover their published internal linear combination 
map shown in Figure 1 (top), one simply applies these 
weights after first multiplying each a\ m by a Gaussian 
beam with FWHM=1° in equation (9). 

The main drawback of this weighting is that it ne- 
glects that there is a tradeoff between foregrounds and 
detector noise which depends strongly on angular scale. 
Diffuse foregrounds are most important on large scales 
where detector noise is negligible, warranting large neg- 
ative and positive weights to aggressively subtract fore- 
grounds. Detector noise, on the other hand, is most im- 
portant on small scales, both because of its Poissonian 
nature (C; roughly constant in the observed map) and 
because the beam correction in equation (10) causes it to 
blow up exponentially on scales smaller than the angular 
resolution [28, 47]. In the limit £ — > oo, the best weight- 
ing is therefore — > (0, 0, 0, 0, 1) regardless of what the 
foregrounds are doing, as illustrated in Figure 4, since 
the W-band has the best resolution. 

We choose to minimize the total unwanted power from 
foregrounds and noise combined, separately for each har- 
monic £ as in [28] , as opposed to only for the combination 
corresponding to the 1° pixel variance as in [3]. As seen 
in Figure 5, one expects such a weighting to combine 
features from the two previous figures: rather aggressive 
subtraction using all five channels at low t, more cautious 
subtraction using only the higher-resolution channels on 
intermediate scales, and all the weight on the W-band 
at extremely high £. In particular, it is crucial to down- 
weight the K-band when cleaning on the scales of the 
acoustic peaks, otherwise the acoustic peaks in the re- 
sulting map will be dominated by K-band noise. 

The constraint equation (5) that we leave the CMB un- 
touched corresponds to the requirement that the weights 
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FIG. 7: Sample band power window functions are shown for the 
cleanest of the sky regions from Figure 8, all normalized so that the 
maximum value is unity. The approximate lack of leakage from odd 
numbers of multipoles away results from the approximate parity 
symmetry of this region. 



where is the 5x5 matrix- valued cross-power spectrum 



= In 1 



(13) 



As an example, Figure 6 shows the weights we obtain for 
the 3rd cleanest of the 9 sky regions shown in Figure 8 be- 
low. We see that just as forecast in Figure 5, and as in the 
WMAP team weighting of Figure 3, the 61 GHz V chan- 
nel is "the breadwinner" , getting a large positive weight 
on large scales since it has the lowest overall foreground 
level. The 94 GHz W channel gets a negative weight to 
subtract out dust, and the three lower frequency channels 
are used to subtract out synchrotron, free-free and any 
other emission dominating at low frequencies. In cleaner 
sky regions, weights get less aggressive in the sense of ac- 
quiring smaller absolute values. In particular, we recover 
weights similar to those of the WMAP team (Figure 3) 
on large scales for the Kp2 sky cut defined and used by 
[3]. 



sum to unity w| = 1) for each £, 



that 



e • w< = 1, 



(11) 



E. Blind analysis and the power spectrum matrix 



where e = (1,1,1,1,1) is a column vector of all ones. 
Minimizing the power (|a^ m | 2 ) in the cleaned map of 
equation (10) subject to this constraint gives [28, 29, 48] 



e*c; 



(12) 



We compute the power spectrum matrix C/ in practice 
using the method of [49] ; a similar approach was used by 
the Boomerang [50] and WMAP [6] teams. This simply 
consists of expanding the masked sky patch in question in 
spherical harmonics and then correcting for window func- 
tion effects. Our only variation is that we do not invert 
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the window matrix to obtain anticorrelated band power 
estimates with delta function window functions. Rather, 
we simply divide by the quadratic estimators by the area 
of the window <5T 2 window functions, which asymptotes 
to the unbiased minimum variance estimators of [51] on 
scales much smaller than the sky patch analyzed. An 
example of our window functions is shown in Figure 7. 

As an example, Figure 9 shows the measured power 
spectrum for the V-band, the cleanest of WMAP's fre- 
quency bands. 

One fact worth emphasizing is that our weighting 
scheme of equation (12) is totally blind, assuming noth- 
ing whatsoever about the CMB power spectrum, the 
foregrounds, the WMAP detector noise or external tem- 
plates. The only assumption is that the CMB spectrum 
is the blackbody that the WMAP team have modeled it 
as, so that the CMB contributes equally to all five chan- 
nels — otherwise the vector e would be replaced by some 
other constant vector. We see that there is no need to 
model the CMB, the foregrounds or the noise, since all 
we need for computing the optimal weights is the total 
power spectrum matrix C;, containing the total contri- 
bution from CMB, foregrounds and noise combined - 
and this can be measured directly from the data. 

We can decompose C; as a sum of two terms, 

C t = Cf nk + Q mb = <T nk + Q mb ee 4 , (14) 

where the second term is the CMB contribution and the 
first term is the contribution from noise and foregrounds. 
Note that if we keep C^ unk fixed and change C^ mh , the 
weights given by equation (12) stay the same. The eas- 
iest way to see this is to note that the quantity we are 
minimizing is (\a£ m \ 2 ) = w'C^w = w t C J £ unk w + C^ mb (e • 
w) 2 = w t C J / ,nk w + Cf mh , so the CMB power is just an 
additive constant that does not affect the optimal weight- 
ing. This means that our method is blind to assumptions 
about the underlying (ensemble-averaged) CMB power 
spectrum. Although we will return below in Section II G 
to the issue of how to determine what fraction of the 
power Ci comes from each of the two terms in equa- 
tion (14), it is important to remember that this affects 
only the physical interpretation, not our cleaning method 
and the maps we produce. 



F. Subdividing the sky 

To minimize the variance in our cleaned map, we 
should take advantage of all ways in which the unwanted 
signals (noise and foregrounds) differ from the CMB in 
their contribution to the covariance matrix N in equa- 
tion (6). Above we exploited their different dependence 
on angular scale t. Unlike the CMB, foregrounds are 
not an isotropic Gaussian random field. Rather, their 
variance differs dramatically between clean and dirty re- 
gions of the sky. It is therefore desirable to subdivide 
the sky into a set of regions of increasing cleanliness 




FIG. 8: The top panel shows our junk map with a color scale 
that is uniform in the logarithm of the temperature in /iK. The 
bottom panel show our subdivision of the sky into seven regions of 
decreasing cleanliness. From outside in, they correspond to junk 
map temperatures T < lOOjUK, 100/LtK - 300ajK, 300/iK - lmK, 
lmK - 3mK, 3mK - lOmK, lOmK - 30mK and T > 30mK, respec- 
tively. The last of these regions contains only the set of roundish 
blobs in the inner Galactic plane. The second dirtiest region is 
seen to be topologically disconnected, and we treat its leftmost and 
rightmost blobs as separate regions, giving nine regions in total. 



and perform the cleaning separately for each one [28]. 
One then expects our method described above to settle 
on more aggressive weights for the dirtier regions, where 
foregrounds are much more of a concern than noise. A 
second advantage of such a subdivision is that the fre- 
quency dependence of the foregrounds is likely to differ 
between very dirty and very clean regions, again result- 
ing in different optimal weights. The WMAP team used 
the latter argument to motivate their subdivision of the 
sky into 12 regions, and convincingly demonstrated that 
foreground spectra indeed did vary across the sky, no- 
tably for synchrotron radiation where the spectrum was 
found to steepen towards increasing galactic latitudes. 

The WMAP team have created a set of sky-masks of 
increasing cleanliness based solely the K-band map. Al- 
though these masks are undoubtedly fine in practice, the 
procedure used in creating them is not blind, since it 
rests on the assumption that all dirty regions are dirty in 
the K-band. In particular, if a foreground manifests it- 
self only at higher frequencies (imagine, say, a blob with 
localized dust emission without detectable synchrotron 
or free-free emission), it would go unnoticed. A second 
minor drawback of this K-band approach is that random 
CMB fluctuations affect the masks at a low level. In 
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other words, the mask was based on the upper left map 
in Figure 2 which, as opposed to the upper right map, 
contains CMB fluctuations. 

To preserve the blind nature of our method, we there- 
fore create sky masks with a different procedure. We first 
form four difference maps W-V, V-Q, Q-K and K-Ka, 
thereby obtaining maps guaranteed to be free of CMB 
signal that pick up any signals with a non-CMB spec- 
trum. We then form a combined "junk map" by taking 
the largest absolute value of these four maps at each point 
in the sky (Figure 8, top). Finally, we create disjoint sky 
regions based on contour plots of this map. We use cuts 
that are roughly equispaced on a logarithmic scale, corre- 
sponding to thresholds of 30000, 10000, 3000, 1000, 300 
and 100/xK (Figure 8, bottom). We emphasize that we 
have found no evidence whatsoever for any actual prob- 
lems with the WMAP team masks, and opt to use our 
own simply to preserve the blind and CMB-independent 
nature of our analysis. As a cross-check, we also repeated 
our entire analysis using the WMAP masks KpO, Lp2 and 
Kpl2, obtaining similar results. 

We followed the WMAP team procedure in the de- 
tails of converting the contours into the masks shown 
in Figure 8 (bottom): we downsampled the junk map to 
HEALPpix resolution 64, imposed the cuts, went back up 
to HEALPix resolution 512, performed Gaussian smooth- 
ing with FWHM=2° on the {0, l}-valued mask and im- 
posed a cutoff of 0.5. The WMAP team reported strong 
spatial variations of foreground spectra in the innermost 
parts of the galactic plane, and therefore subdivided this 
into 11 disjoint regions. We were unable to reproduce 
this procedure since they did not specify which these re- 
gions were. Instead, we merely lopped off three spatially 
disconnected islands in the two dirtiest regions as their 
own separate masks, as illustrated in Figure 8 (bottom), 
leaving 9 separate masks in total. 

Our multipole-based cleaning is nonlocal in the spirit of 
equation (4), and although it guarantees that the CMB 
signal is preserved separately for each pixel, this is of 
course not the case for foregrounds. To avoid mirages of 
foreground emission from the Galactic plane leaking up 
to high latitudes, we clean the galaxy "from inside out" , 
i.e., clean the dirtiest region first, the second dirtiest re- 
gion second, etc. To be specific, we start by defining five 
temporary maps, initially set equal to the five WMAP 
channels (Figure 2, left). We then repeat the follow- 
ing cleaning procedure nine times, once corresponding to 
each region i = 1, ...,9: 

1. Compute the power spectrum matrix Ci and the 
optimal weights for the i th region only. 

2. Expand the five temporary maps in spherical har- 
monics and compute a cleaned all-sky map using 
the weights from step 1. 

3. Replace the i th region of the temporary maps by 
the corresponding region in the cleaned map from 
step 2 (smoothed to the resolution appropriate for 



that channel, of course, by multiplying by B\ in 
^-space). 

Visually, one thus sees the foreground contamination 
gradually being cleaned off from Figure 2 as this iter- 
ation proceeds, starting in the inner Galactic plane and 
proceeding outward. At the end of the nine iterations, 
the five temporary maps equal the desired cleaned map 
at the K-, Ka-, Q-, V-, and W-band resolutions, respec- 
tively. 



G. Interpretation of the cleaned maps 

Figure 1 (bottom) shows our final cleaned map and 
Figure 10 shows its power spectrum in sky regions of 
varying cleanliness. We plot all maps after 21' Gaussian 
smoothing (giving a net FWHM of 24') to prevent them 
from being undersampled by the pixels in the image. The 
reader interested in using this map can download the 
corresponding 13 Megabyte HEALPix fits file from the 
web 3 . To use this map, it is important to be clear on 
how to interpret it. 

First of all, it is a sum of CMB, foreground and de- 
tector noise fluctuations. Although it was constructed 
by minimizing its power spectrum Cf can = w^C^w^, 
its power nonetheless gives a lower bound on the CMB 
power spectrum, since the minimization was performed 
subject to the constraint that the CMB be preserved; 
e • w = l 4 . 

So what fraction of the power seen in Figure 10 is due 
to CMB, foregrounds and detector noise, respectively? 
Let us first get some rough estimates from the figures, 
then present more quantitative limits. Since this paper 
is focused on minimizing foregrounds, not on physically 
modeling them, the interested reader is referred to the de- 
tailed foreground study of the WMAP team [3] for more 
information. 

The WMAP team cleverly eliminated the average noise 
contribution by using only cross-correlations between 
different differencing assemblies (DAs) to measure the 
power spectrum, and used a combination of sky cutting 
and foreground subtraction (with external templates) to 
minimize the foreground contribution. The best fit cos- 
mological model [17, 18] to their measured CMB power 
spectrum [6] is shown for comparison in figures 9 and 10. 



3 Our cleaned maps are available for download at http://www. 
hep.upen.edu/-max/wmap.html, together with a high-resolution 
version of this paper. 

4 The only way in which its power could be biased low would be if 
random fluctuations in our estimate of the Ci matrix conspired 
to remove power. Although we found no indication of this ac- 
tually happening, we computed our weights w< using a heavily 
smoothed version of C^-matrix as a precaution. Specifically, we 
smooth over at least At = 10 or 100 (£, m)-modes, whichever is 
larger, obtaining an ^-dependence of Ci with no visible trace of 
random fluctuations. 
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Multipole 1 

FIG. 9: Total power spectra for CMB, foregrounds and noise com- 
bined are shown for the 61 GHz V-band, the one with the overall 
lowest foreground levels. This is the (4, 4) element of the power 
spectrum matrix C^. From bottom to top, they correspond to the 
five cleanest sky regions shown in Figure 8. For comparison, the 
thick curve shows the best-fit CMB model from [17, 18]. 




Multipole 1 

FIG. 10: Same as previous figure, but for our foreground-cleaned 
CMB map (bottom panel of Figure 1). The power spectra of the 
cleanest sky regions are seen to be virtually identical with those 
for V-band on large scales, showing how subdominant smooth fore- 
grounds arc at 61 GHz. 
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Multipole 1 



FIG. 11: Power spectra of non-CMB signals (foregrounds plus de- 
tector noise) in the cleanest part of the sky shown in Figure 8. The 
five thin curves show the power spectra of five maps on the right 
side of Figure 2, i.e., the five WMAP channels after our cleaned 
CMB map (Figure 1, bottom) has been subtracted out. From top 
to bottom, the five curves correspond to 22.8, 33.0, 40.7, 93.5 and 
60.8 GHz, respectively (note that the second highest frequency, 
V-band, is the cleanest). These curves should be interpreted as 
lower limits on foregrounds plus noise. The black curve gives the 
lower limit on foregrounds plus noise in our cleaned map using the 
method described in the text. 



For I < 100, it is seen to agree well with the lower en- 
velope of our curves in Figure 9, suggesting that fore- 
grounds are subdominant in the cleanest parts of the 61 
GHz sky. Figure 10 shows no noticeable excess power 
at U < 100 due to foregrounds in the cleaned map in 
any of the four cleanest sky regions from Figure 8, which 
together cover all but the very innermost Galactic plane. 

The slight power deficit on the very largest scales has 
two causes: one is the low quadrupole, to which we return 
in Section III below, which pulls down neighboring band 
power estimates because the band-power window, illus- 
trated in Figure 7, contains small off-diagonal contribu- 
tions. The second cause is that we have not corrected for 
the effects of monopole and dipole removal, which pulls 
down the power estimates on the scale of the sky patch in 
question (our method produces an unbiased CMB map 
regardless of what Cg we use in equation (12), so this 
merely raises the variance in the cleaned map above the 
optimal level). On smaller scales, detector noise starts 
to dominate, and is seen to push our curves way above 
the CMB curve. A worthwhile future project for further 
quantifying foregrounds would be to repeat our analysis 
with C( estimated in a way that removes the detector 
noise contribution. For Q, V and W bands, which each 
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have more than one DA, this can be done using cross- 
correlations. For the (1, 1) and (2, 2) components of 
(the K and Ka power spectra), it would involve subtract- 
ing the noise power using the WMAP team's noise model. 

We can, however, give some quantitative limits even 
based on our measured Ci alone. Grouping the five co- 
efficients a\ m into a five-dimensional vector a^ m , equa- 
tion (13) becomes simply — (a| TO a^ m ), and the 
cleaned procedure can be written ae m = w • a£ m . Us- 
ing equation (12) shows that the power in the cleaned 
map is 

CV cloan = (15) 

By subtracting our cleaned map from the input map, 
we obtain maps showing in Figure 2 (right). These are 
guaranteed to be free from CMB power, since the cleaned 
map gave a lower limit on the CMB. Most likely, we have 
subtracted some foreground power too, so the five maps 
should be interpreted as placing lower limits on the fore- 
ground power. (As a toy example, imagine synchrotron, 
free-free and dust emission tracing each other perfectly; 
the sum of their three spectra can then be written as a 
constant, which our method will interpret as CMB, plus a 
non-negative residual, which our method will interpret as 
foregrounds). Using equation (14), the covariance matrix 
of these five CMB-free "junk maps" is 

C J / nk > Q ocmb = C t - C| lean ee* = II<C<II|, (16) 

where the projection matrix 

U f = I ewj (17) 

satisfies nf = Tig, Tlge = and II^w = and can be 
interpreted as projecting out the CMB component. In 
this notation, the maps in Figure 2 (right) are defined 
by simply U e a em . The inequality Cf nk > C£ ocmb refers 
only to the diagonal elements of these matrices. The five 
diagonal elements of C^ unk are plotted in Figure 11 for 
our cleanest sky region from Figure 8. 

How much of can possibly be due to CMB? In other 
words, forgetting for a moment about the weighting that 
gave equation (15), how large can we make C| mb in equa- 
tion (16) before the covariance matrix C J ^ un gets unphys- 
ical properties? First of all, no covariance matrix can 
have negative eigenvalues, so we we must stop increasing 
C£ mb once the smallest eigenvalue of C J £ unk drops to zero. 
In fact, n*w = implies that when using our optimal 
weighting and hence equation (15), C^ unk has a vanish- 
ing eigenvalue corresponding to the vector w, and it is 
easy to show that this alternative method for estimating 
C£ mb is equivalent to our original method, being simply 
an alternative derivation of equation (15). 

Let us now make a second assumption: that C^ unk 
cannot have any negative elements. The noise covari- 
ance matrix is guaranteed to have this property, since 
the absence of correlations between bands implies that 



it is diagonal. The foreground covariance matrix is also 
guaranteed to have this property if foreground emission is 
indeed emission, i.e., if foregrounds can make only posi- 
tive contributions to the sky maps. Pure absorption at all 
frequencies likewise give only positive correlations. This 
assumption, made also in the Maximum-Entropy analy- 
sis of [3] is likely to be valid WMAP frequencies, since 
the only known exception is the thermal SZ effect, and it 
changes sign only outside of the WMAP frequency range 
(around 217 GHz). In summary, we therefore have two 
separate limits on the CMB power spectrum: 

C| mb < Cf e&n , (18) 
Q cmb < nrin(Q)y. (19) 

£iciean j g a j g0 ^ e ac t ua l variance of the cleaned map, so 
if the second limit is lower than first, then we know that 
the difference cannot be due to CMB. 

To gain intuition, consider the following two examples 
of the junk + CMB decomposition of equation (14) for 
the simple case of only two frequency bands: 




In the first case, equation (15) gives C| mb = 1, i.e., 
the largest contribution we can possible attribute to the 
CMB (the second term) is ee* — if we scaled up the 
last term, then C J ^ unk (the first term) would acquire an 
unphysical negative eigenvalue. Moreover, we see that 
C^ unk looks like the covariance matrix of a perfectly re- 
spectable foreground, with a spectrum such that it con- 
tributes twice as high rms fluctuations to the first channel 
as to the second, and with perfect correlation between the 
two (the dimcnsionless correlation coefficient is r = 1). 
Such a foreground would be completely removed by our 
method, and so we cannot place any lower limit on the 
foreground contribution to the cleaned map in this case. 
In the second case, equation (15) gives C£ mb = 3/2, 
so this will be the variance in the cleaned map. How- 
ever, C^ unk (the first term) in equation (21) exhibits an 
unphysical anticorrelation between the two bands which 
neither noise nor foregrounds could produce, and equa- 
tion (19) tells us that the CMB power C| mb < 1, cor- 
responding to the alternative decomposition on the last 
line. This means that the junk map must have a power 
contribution of at least 1/2 which is not due to CMB. 

This difference, which places a lower limit on the 
amount of non-CMB signal in our cleaned map, is shown 
as the black line in Figure 1 1 for our actual five-band case. 
We see that although we get an interesting lower limit for 
I > 100, presumably dominated by detector noise, the 
residual foreground level is consistent with zero at low I. 
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H. Wiener filtering 

Figure 1 (bottom) shows our cleaned map of the CMB 
sky, and this is the map that should be used for cross- 
correlation analysis with other data sets and other sci- 
entific applications. For visualization purposes, however, 
we can do better. The "best-guess" map of what the 
CMB looks like, in the sense of minimizing the rms er- 
rors, is the Wiener filtered map defined by [52] 



Wiener 



c: 



cmb 



C: 



clean 



0>£n 



(23) 



This staple signal processing technique, multiplying by 
"signal over signal plus noise" , has additional attrac- 
tive properties; for instance, it constitutes the best-guess 
(maximum posterior probability) map in the approxima- 
tion of Gaussian fluctuations. Examples of recent ap- 
plications of Wiener filtering to CMB and galaxy map- 
ping include [53-59]. Our resulting Wiener filtered map 
is shown in Figure 1 (middle), using the best fit model 
from the WMAP team [17, 18] as our estimate of C| mb 
in the numerator of equation (23). For the denomina- 
tor, we take the larger of our measured C| lean and C| mb , 
so that the ratio is guaranteed to be < 1. The result 
is not an unbiased CMB map. Rather, equation (23) 
shows that each multipole gets multiplied by a number 
between and 1, so features with high signal-to-noise are 
left unaffected whereas features that are not statistically 
significant become suppressed. This means that the fea- 
tures that you see in Figure 1 (middle) are likely to be 
real CMB fluctuations, having signal-to-noise exceeding 
unity. 

The cleaned map at the bottom of Figure 1 reveals 
some residual galactic fluctuations on very small angular 
scales, caused mainly by the fact that no other channels 
have fine enough angular resolution to help clean the W- 
band map for very large I. We Wiener filter each of 
the nine sky regions from Figure 8 separately, so Cf e&n 
is much higher near the Galactic plane. This is why the 
Galactic contamination is imperceptible in the Wiener fil- 
tered map: equation (23) automatically suppresses fluc- 
tuations in regions with large residual foregrounds. 



III. DISCUSSION 

We have performed an independent foreground analy- 
sis of the WMAP maps to produce a cleaned CMB map. 
The only assumption underlying our method is that the 
CMB contributes equally to all five channels. This as- 
sumption rests on very solid ground [1] . The basic reason 
for this is that the COBE/FIRAS determination of the 
CMB spectrum [60] is based on the absolute CMB sig- 
nal, which is about 10 5 times larger than the fluctuations 
that we have considered in this paper. 

Figure 1 shows that our map agrees very well with the 
internal linear combination (ILC) map from the WMAP 




FIG. 12: Comparison of the total (CMB+foregrounds+noise) 
power spectra of the WMAP team internal linear combination map 
[3] (top curve) and our cleaned map (middle curve), both for our 
cleanest sky region. Both of these cleaned maps are seen to re- 
produce on large scales the CMB power spectrum measured by 
the WMAP team [6] (lower wiggly curve), which has no net noise 
contribution because it is based on cross-correlations between chan- 
nels. The lower smooth curve is the WMAP team's best fit model 
[f7, 18]. As explained in the text, the noise contribution is seen 
to become important earlier in the WMAP team's cleaned map 
than in ours because it is limited by the lowest resolution frequency 
bands. All power spectra have been smoothed with a boxchar filter 
of width A£ = 10 to reduce scatter. 



team on the scales > 1° where they can be compared. 
This is yet another testimony to the high signal-to-noise 
in the WMAP data and to the fact that unpolarized CMB 
foregrounds are manageable: the basic spatial features of 
the CMB are insensitive to the details of the foreground 
removal method used. 

The basic advantage of our map is illustrated in Fig- 
ure 12, which compares its all-sky power spectrum with 
that of the WMAP team ILC map. Both power spec- 
tra have had beam effects removed here; in Figure 1 , our 
map is shown at the W-band resolution and the ILC map 
is shown at 1° resolution as released. 

The first thing to note from Figure 12 is that fore- 
grounds appear highly subdominant in both maps, since 
their power spectra essentially coincide with the WMAP 
CMB power spectrum on large scales where noise be- 
comes unimportant. Second, as expected, the main im- 
provement in our map is seen to be on the smaller scales 
where noise is important, gaining a factor of 30% at the 
first acoustic peak and about a factor of two at the sec- 
ond peak where noise from the low frequency channels is 
beginning to exponentially dominate the ILC map. 
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FIG. 13: A heretic all-sky analysis of our cleaned map (solid jagged 
curve) gives a slightly less low quadrupole <5T 2 than the cut-sky 
WMAP analysis [6] (dashed curve), and also agrees well with the 
quadrupole from an all-sky analysis of the WMAP team cleaned 
map of [3] (dotted curve). The smooth curve shows the WMAP 
team best fit model [17, 18] with the band indicating the cosmic 
variance errors (WMAP noise and beam effects are completely neg- 
ligible on these scales). 



FIG. 14: The angular correlation functions are shown corre- 
sponding to our all-sky cleaned map (solid curve), the cut-sky 
WMAP power spectrum (dashed curve) and the best-fit WMAP 
team cosmological model (dotted curve). This is simply the Leg- 
endre transform of the power spectra from the previous figure, i.e., 
(47r) _1 Y2t @£ + 1)-Pf (cos 9)Ci, so the the quadrupole and octopolc 
contributions are shaped as 3 cos 2 6 — 1 and 5 cos 3 9 — 3 cos a;, re- 
spectively. 



We hope that our map will prove useful for a vari- 
ety of scientific applications. For cross-correlation with 
external maps, its lowered noise power should be partic- 
ularly advantageous for pulling out small-scale signals, 
for instance those associated with lensing and SZ clus- 
ters. Rather that attempting detailed modeling of the 
residual noise and foreground fluctuations in our map, a 
simple way to place error bars on such correlations will be 
repeating the analysis with a suite of rotated and flipped 
versions of our map as in, e.g., [61]. 

Table 1 - Measurements of the CMB quadrupole and octopolc. 



Measurement 


<5T 2 [ M K 2 ] 


p- value 


«5T| [/iK 2 ] 


Spergel et al. model 


869.7 




855.6 


Hinshaw et al. cut sky 


123.4 


1.8% 


611.8 


1LC map all sky 


195.1 


4.8% 


1053.4 


Cleaned map all sky 


201.6 


5.1% 


866.1 


Cosmic quadrupole 


194.2 


4.7% 




Dynamic quadrupole 


3.6 







Let us close by returning from small to large angu- 
lar scales. The surprisingly small CMB quadrupole has 
intrigued the cosmology community ever since if was 
first observed by COBE/DMR [62], and simulations by 
[17, 18] show that the low value observed by WMAP 
is sufficiently unlikely to warrant serious concern. The 



WMAP team measured the quadrupole using only the 
part of the sky outside of their Galactic cut, and stress 
that the dominant uncertainty in its value is foreground 
modeling. While we agree with this assessment, it should 
be borne in mind that noise variance and beam issues 
are completely negligible on these huge angular scales, so 
this does not imply that the foreground uncertainties are 
large compared to the signal itself. Indeed, Figure 12 sug- 
gests that foregrounds are subdominant to the intrinsic 
CMB signal even without any Galaxy cut if a foreground- 
cleaned map is used. While the reader may feel disturbed 
by the clearly visible Galactic residuals in Figure 1 (in 
both the top and bottom maps), it is important to bear 
in mind that these signals are present in only a tiny frac- 
tion of the total sky area and therefore contribute little 
to the total power spectrum. 

To quantify the impact of the Galactic plane, we com- 
puted a sequence of all-sky power spectra for both our 
map and the ILC map. We first used the unadulterated 
maps, then repeated the calculation after repacing the 
dirtiest of the nine regions from Figure 8 with zeroes, 
then zeroed our the 2nd dirtiest region as well, etc. This 
zeroing procedure obviously biases the measurements by 
removing power, by an amount related to the zeroed 
area, but provides a powerful test of how sensitive the 
results are to Galactic plane details. We also band-pass 
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FIG. 15: The left panel shows the quadrupole (top), octopole (middle) and hexadecapole (bottom) components of our cleaned all-sky 
CMB map from Figure 1 on a common temperature scale. Note that not only is the quadrupole power low, but both it and the octopole 
have almost all their power perpendicular to a common axis in space, as if some process has suppressed large-scale power in the direction 
of this axis. We computed the corresponding images for the WMAP team ILC map as well, and found them to be very similar. The 
right panel shows he cosmic quadrupole (top) after a correcting for a crude estimate of the dynamic quadrupole (middle) from our motion 
relative to the CMB rest frame. The bottom right map shows the sum of the quadrupole and octopole maps from the left panel. 



filtered the resulting maps to produce spatial plots of the 
quadrupole, octopole and hexadacapole. We found that 
zeroing out the dirtiest parts of the Galactic plane had a 
negligible effect on both the power spectrum and on the 
spatial structure of these lowest multipoles. Specifically, 
we could zero out all but the three cleanest regions Fig- 
ure 8 (everything with T > lmK) without the quadrupole 
or octopole changing substantially. The spatial morphol- 
ogy of the quadrupole, octopole and hexadecapole for the 
all-sky analysis agrees well between our map and the ILC 
map, again showing insensitivity to galaxy modeling de- 
tails (in particular, the ILC map is likely to have less 
contamination in the Galactic plane due to more subdi- 
visions there). 

Although more detailed foreground modeling would be 
needed to rigorously quantify the foreground contribu- 



tion to low multipoles, let us, encouraged by the above- 
mentioned tests, tentatively assume that this contribu- 
tion is unimportant and perform an all-sky analysis of 
our cleaned map. The resulting power spectrum is shown 
in Figure 13, which is simply a blow-up of the leftmost 
part of the previous figure, and the corresponding angu- 
lar correlation function is shown in Figure 14. 

Table 1 summarizes the quadrupole and octopole re- 
sults. We see that although the quadrupole is still low, 
it is not quite as low as that from the cut-sky WMAP 
team analysis of [6] . Moreover, our map has a quadrupole 
virtually identical to the WMAP team ILC map despite 
the differences in foreground modeling, further support- 
ing our hypothesis that the quadrupole is not strongly 
affected by foregrounds. The second column in Table 1 
shows the probability of the quadrupole in our Hubble 
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volume being as low as observed if the best-fit WMAP 
team model from [17, 18] is correct. This is computed 
for the all-sky case where 5T^ has a ^-distribution with 
5 degrees of freedom, so the probability tabulated is sim- 
ply 1 -7[5/2, (5/2)T 2 2 /855.6/xK 2 ]/r(5/2), where 7 and V 
are the incomplete and complete Gamma functions, re- 
spectively 5 . We see that the statistical significance of the 
low quadrupole problem drops substantially with our all- 
sky analysis, below the 95% significance level, requiring 
merely a one-in-twenty fluke. 

To understand why the all-sky quadrupole is larger 
than the cut-sky one measured by the WMAP team, we 
plot the lowest three multipoles of our cleaned map in 
Figure 15, all on the same temperature scale. Several 
features are noteworthy. First of all, the quadrupole is 
low with an interesting alignment. A generic quadrupole 
has three orthogonal pairs of extrema (two maxima, two 
minima and two saddle points). We see that the actual 
CMB quadrupole has its strongest pair of lobes, appar- 
ently coincidentally, fall near the Galactic plane. Apply- 
ing a Galaxy cut therefore removes a substantial fraction 
of the quadrupole power. The saddle point is seen to be 
close to zero. In other words, there is a preferred axis in 
space along which the observed quadrupole has almost 
no power. 

Second, the observed quadrupole is the sum of the cos- 
mic quadrupole and the dynamic quadrupole due to our 
motion relative to the CMB rest frame - see [63] for a de- 
tailed discussion. Since this motion is accurately known 
from the CMB dipole measurement [1, 62], the dynamic 
quadrupole can and should be subtracted when studying 
the cosmic contribution. Figure 15 shows the dynamic 
quadrupole approximated by (v/c) 2 (cos6 2 — 1/3), where 
v s=a 369km/ s towards (l,b) « (264,48) is the velocity 
of the Solar System relative to the CMB [1] and 6 is 
the angle relative to this velocity vector. This is a small 
correction with peak-to-peak amplitude {v/c) 2 ~ 4^(K, 
and Table 1 shows that it reduces the cosmic quadrupole 
slightly. This approximation is crude since the dynamic 
quadrupole in fact depends on frequency [63], and there- 
fore may have been cither over- or under-estimated in our 
foreground cleaned map - we have shown it here merely 
to illustrate its spatial orientation and give an order-of- 
magnitude estimate of its importance. 

Third, although the overall octopole power is large, 
not suppressed like the quadrupole, it too displays the 
unusual property of a preferred axis along which power 
is suppressed. Moreover, this axis is seen to be approx- 
imately aligned with that for the quadrupole. The rea- 
son that our measured octopole in Figure 13 is larger 
than that reported by the WMAP team is therefore, once 



5 The actual probability is slightly larger for the cut-sky case where 
there are fewer effective degrees of freedom, although [17] show 
with Monte Carlo simulations that there it is still a disturbingly 
small for the measured WMAP quadrupole — indeed, they ob- 
tain a consistency probability of 0.7% using Markov chains. 



again, that much of the power falls within the Galaxy cut. 
In contrast, the hexadecapole is seen to exhibit the more 
generic behavior we expect of an isotropic random field, 
with no obvious preferred axis. 

How significant is this quadrupole-octopole alignment? 
As a simple definition of preferred axis for an arbi- 
trary sky map, [64] computes the unit vector r around 
which the angular momentum dispersion ((? • L) 2 ) = 
m2 \ a im(j)\ 2 is maximized. Here a^ m (r) denotes the 
spherical harmonic coefficients of the map in a rotated 
coordinate system with its z-axis in the the r-direction. 
The preferred axes f 2 and f 3 for the quadrupole and oc- 
topole, respectively, are 

? 2 = (-0.1145,-0.5265,0.8424), 

r 3 = (-0.2578,-0.4207,0.8698), (24) 

i.e., both roughly in the direction of (I, b) ~ (—110°, 60°) 
in Virgo. The angle between these two axes is merely 
10.3°, and the probability that a random axis falls inside 
a circle of radius 10.3° around the quadrupole axis is 
simply the area of this circle over the area of the half- 
sphere, about 1/62. In other words, if the CMB is an 
isotropic Gaussian random field, then a chance alignment 
this good requires about a l-in-62 fluke. This issue is 
discussed in greater detail in [64] . 

What does this all mean? Although we have presented 
these low multipolc results merely in an exploratory 
spirit, and more thorough modeling of the foreground 
contribution to I = 2 and I — 3 is certainly warranted, 
it is difficult not to be intrigued by the similarities of 
Figure 13 with what is expected in some non-standard 
models, for instance ones involving a flat "small Uni- 
verse" with a compact topology as in [65-70] and one 
of the three dimensions being relatively small (of order 
the Horizon size or smaller). This could have the effect of 
suppressing the large-scale power in this particular spa- 
tial direction in the same sense as is seen in Figure 13. 
As so often in science when measurements are improved, 
WMAP has answered old questions and raised new ones. 
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